/*
#========================================================================
#
# Alan Rogers
# Massachusetts Institute of Technology Haystack Observatory
# Program: vlbi2
# Calculate cross or autocorrelation spectra from Mk5b recorded data.
#
# The vlbi2 program is used to correlate two Mark 5B files with the same
# starting record time and record mode. A postscript file dd1.pos is
# produced showing the correlation amplitude and phase as a function
# of frequency in all channels in the two files with all channels
# plotted on the same axes. Correlation coefficients (fringe amplitudes,
# which will be 1.00 for an autocorrelation) are indicated below each
# channel. The vlbi2 program accepts two files, each 10 MB in size,
# and two flags. The usage is:
#
#        vlbi2 <file1> <file2> [-2bit <2bit_arg>] [-rev <rev_arg>]
#
# Where the items in [] are optional. The -2bit flag indicates whether
# the files are 2 bits/sample recordings (2bit_arg is 1) or 1 bit/sample
# (2bit_arg is 0) [default]. The files must be the same number of
# bits per sample. The .rev flag determines whether the channels in dd1.pos
# are plotted in reverse order (rev_arg is 1) or normal order
# (rev_arg is 0) [default]. It is important to note that with the .rev
# flag, though the channels can be plotted in the reverse order, they
# are not flipped in frequency. If file1 and file2 are the same, an
# autocorrelation is produced; if they are different a cross correlation
# is produced.
#
#
#=======================================================================
#
# Dec 2010
# S.Weston              stuart.weston@aut.ac.nz
# AUT University
#
# Modified to use DISLIN so that we can see the data better than just
# printed across the bottom of an A4 sheet in postscript. Also use
# DISLIN widgets to provide a better UI.
#
#=======================================================================
*/

#include <stdio.h>
#include <math.h>
#include "vlbi2.h"

void
Four (double *fft_r, double *fft_i, int nn)
{
  int n, mmax, m, j, istep, i, ii, jj;
  double wtemp, wr, wpr, wpi, wi, theta;
  double tempr, tempi;

  n = nn << 1;
  j = 1;
  for (i = 1; i < n; i += 2)
    {
      if (j > i)
	{
	  jj = (j - 1) / 2;
	  ii = (i - 1) / 2;
	  tempr = fft_r[jj];
	  tempi = fft_i[jj];
	  fft_r[jj] = fft_r[ii];
	  fft_i[jj] = fft_i[ii];
	  fft_r[ii] = tempr;
	  fft_i[ii] = tempi;
	}
      m = n >> 1;
      while (m >= 2 && j > m)
	{
	  j -= m;
	  m >>= 1;
	}
      j += m;
    }
  mmax = 2;
  while (n > mmax)
    {
      istep = mmax << 1;

      theta = -(6.28318530717959 / mmax);
      wtemp = sin (0.5 * theta);
      wpr = -2.0 * wtemp * wtemp;
      wpi = sin (theta);
      wr = 1.0;
      wi = 0.0;
      for (m = 1; m < mmax; m += 2)
	{
	  for (i = m; i <= n; i += istep)
	    {
	      j = i + mmax;
	      jj = (j - 1) / 2;
	      ii = (i - 1) / 2;
	      tempr = wr * fft_r[jj] - wi * fft_i[jj];
	      tempi = wr * fft_i[jj] + wi * fft_r[jj];
	      fft_r[jj] = fft_r[ii] - tempr;
	      fft_i[jj] = fft_i[ii] - tempi;
	      fft_r[ii] += tempr;
	      fft_i[ii] += tempi;
	    }
	  wr = (wtemp = wr) * wpr - wi * wpi + wr;
	  wi = wi * wpr + wtemp * wpi + wi;
	}
      mmax = istep;
    }

}
